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Abstract 



The adaptation of large asexual populations is hampered by the competition between in- 
dependently arising beneficial mutations in different individuals, which is known as clonal 
interference. In classic work, Fisher and MuUer proposed that recombination provides an 
evolutionary advantage in large populations by alleviating this competition. Based on recent 
progress in quantifying the speed of adaptation in asexual populations undergoing clonal in- 
terference, we present a detailed analysis of the Fisher-Muller mechanism for a model genome 
consisting of two loci with an infinite number of beneficial alleles each and multiplicative 
(non-epistatic) fitness effects. We solve the deterministic, infinite population dynamics ex- 
actly and show that, for a particular, natural mutation scheme, the speed of adaptation 
in sexuals is twice as large as in asexuals. This result is argued to hold for any nonzero 
value of the rate of recombination. Guided by the infinite population result and by previous 
work on asexual adaptation, we postulate an expression for the speed of adaptation in finite 
sexual populations that agrees with numerical simulations over a wide range of population 
sizes and recombination rates. The ratio of the sexual to asexual adaptation speed is a 
function of population size that increases in the clonal interference regime and approaches 
2 for extremely large populations. The simulations also show that recombination leads to a 
strong equalization of the number of fixed mutations in the two loci. The generalization of 
the model to an arbitrary number L of loci is briefly discussed. For a particular communal 
recombination scheme, the ratio of the sexual to asexual adaptation speed is approximately 
equal to L in large populations. 
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The evolutionary advantage of sex remains one of the most intriguing puzzles in evolu- 
tionary biology ( IKoNDRASHOV 19931 Ide Visser and Elena 2007| IQtto 2009j) . Many hy- 



potheses have been suggested explaining why sexual reproduction is widespread in nature de- 
spite apparent disadvantages such as the two- fold cost of sex (IMaynard Smith 1978p . Well- 
known examples are the deterministic mutation hypothesis (IKondrashov 1988p . the Fisher- 
MuUer mechanism (IFlSHER 1930[ IMuller 1932[ ICrow and KiMURA 1965P and MuUer's 



ratchet (IMuller 1964[ IFelsenstein 1974|) . to name only a few. These three hypotheses 
are applicable when the fitness landscape in question has certain specific features. Specifi- 
cally, the deterministic mutation hypothesis requires deleterious mutations to be synergis- 
tically epistatic, while the Fisher-MuUer (FM) mechanism as well as MuUer's ratchet can 
explain the advantage of sex if epistasis is negligible. 

Theoretical analyses of the effect of epistasis on the speed of MuUer's ratchet have con- 
cluded that it practically stops operating when epistasis is synergistic (ICharlesworth et al. 1993t 
IKondrashov 1994[IJain 2008p . Furthermore, recent experimental analyses of empirical fit- 



ness landscapes seem to indicate that a particularly strong form of epistasis termed sign epis- 
tasis ( IWeinreich et al. 2005P is the rule rather than an exception (IWeinreich et al. 2006| 



Ide Visser et al. 2009| IFranke et al. 201l| ISzendro et al. 2013|) . Sign epistasis gener- 



ally implies that the fitness landscape is rugged. On a rugged fitness landscape sex can 
be detrimental, even without taking into account the two-fold cost of sex, in that sexual 
populations, unlike the corresponding asexual populations, cannot escape from local fit- 
ness peaks (ICrow and Kimura 1965| IEshel and Feldman 1970| Ide Visser et al. 2009] 
IPark and Krug 20111) . 



Although research on empirical fitness landscapes has been growing substantially in recent 
years, it is still practically infeasible to reliably determine genotypic fitness on a genome- 



wide scale [but see KoUYOS et al. (2012)]. Because of the small sizes of most empirical 
fitness landscapes that have so far been constructed experimentally, the implications of 
the ubiquity of sign epistasis observed in these studies for long term evolution remain un- 



clear. At the same time experimental evidence in favor of the FM mechanism has accumu- 
lated fICOLEGRAVE 2002t ICOOPER 20071) . For these reasons further quantitative analysis of 



the advantage of sex in the absence of epistasis remains a worthwhile endeavor, and we will 
pursue this approach in the present contribution. 

The essence of the FM mechanism is the competition between independently arising bene- 
ficial mutations, termed clonal interference, which slows down the adaptation of large asexual 
populations (IGerrish and Lenski 1998t IMiralles et al. 1999| IWilke 20"04l IKiM and Orr 20051 
IPark and Krug 20071 IFogle et al. 20081ISniegqwski and Gerrish 20T0tlScHiFFELS et al. 201 ip . 



The concept of clonal interference has played an important role in intepreting the behavior 
observed in laboratory selection experiments (ILenski et al. 199lt ILenski and Travisano 1994} 
IBarrick et al. 2009p . and has also been invoked in explaining the population-size depen- 



dence of evolutionary predictability in rugged fitness landscapes (IJAIN et al. 201H ISzendro et al. 2013P 



Although in its original formulation clonal interference theory neglects the occurrence of sec- 
ondary beneficial mutations within a growing clone (IGerrish and Lenski 1998| IGerrish 200ip 



in general the coexistence of multiple beneficial mutations cannot be neglected in large pop- 
ulations (IPark and Krug 2007p . In the following we will therefore use the term clonal 
interference in a wider sense than originally conceived, in that two clones with different 
numbers of beneficial mutations can compete with each other for fixation. 

Much recent theoretical work has focused on obtaining accurate quantitative estimates 
of the speed of adaptation in the presence of clonal interference for the simple situation 
of an unlimited supply of beneficial mutations that act independently on fitness, without 



epistatic interactions [see Park et al. (2010) for review]. It turns out that the population 
dynamics in this regime is well described by a traveling wave moving at constant speed along 
a one-dimensional fitness space. The traveling wave picture was first established for the case 



when the beneficial selection coefficient is the same for all mutations (ITsimring et al. 1996| 



IRouziNE et al. 2003| IDesai and Fisher 20071 IBrunet et al. 20081 IRouzine et al. 2008p 



and recently extended to the more realistic case of selection coefficients drawn from a con- 
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tinuous effect size distribution flGoOD et al. 2012} TFisher 20131) . 

In natural populations it is unlikely that the traveling wave picture persists forever. Apart 
from the assumed absence of epistatic interactions, there are two main features that lead to a 
breakdown of this picture in long-term evolution. First, a fluctuating environment generally 
makes the fitness landscape change with time. In a time-dependent situation it is problematic 
to compare fitnesses of two individuals living on different landscapes and the notion of 



absolute fitness has to be replaced by fitness flux (IMUSTONEN and Lassig 2010j) . Second, 
even if the fitness landscape remains constant for a very long time, the indefinite supply 
of beneficial mutations appearing at constant rate cannot be a good approximation in the 
real world. For example, in long-term evolution experiments the speed of adaptation usually 
slows down (ILenski and Travisano 19941 IBarrick et al. 200"9|) . which is attributed to 



the decreasing supply of beneficial mutations. In this context, the house-of-cards model, 
in which fitness values are assigned randomly to genotypes, could provide a more realistic 
description ( IKingman 1978t IPark and Krug 2008^ . In the framework of this model one 
cannot however explain the advantage of sex, because recombination is no different from 
mutation. 

Although the non-epistatic model with an infinite supply of beneficial mutations is of 
limited validity, it can provide a reasonable approximation when a population undergoes a 
severe environmental change, as is often the case at the beginning of an evolution experiment. 
At the same time this setting is conceptually simple and allows for detailed (if approximate) 
mathematical analysis. In the present paper, we therefore build upon the recent line of work 
on asexual populations undergoing clonal interference and add to it a minimal yet realistic 
recombination scheme. 

Specifically, we consider a facultatively sexual population model with two genetic loci, 
each of which can acquire infinitely many beneficial mutations. We assume that epistasis 
is absent both between and within loci. We first analyze the infinite population dynamics, 
obtaining exact expressions for the speed of adaptation in the limiting cases of zero and 

6 



maximal recombination rates (asexuals vs. obligate sexuals). When the selection coefficient 
of beneficial mutations is the same at both loci and at most one mutation may occur per 
generation and individual, the speed of adaptation for obligate sexuals is twice that of 
asexuals, a result that we argue holds for any positive recombination rate. Based on this 
observation and on previous results for asexual populations, we propose a simple expression 
for the speed of adaptation in finite sexual populations that is in good agreement with 
simulations and predicts a roughly 'two-fold advantage of sex' for large populations and 
moderate recombination rates. In DISCUSSION we describe a possible extension of the 
model to more than two loci, and place our work into the context of related studies. 

MODELS 

We consider a (sexual or asexual) population of haploid individuals in discrete generations. 
The population size is denoted by N and assumed to be constant. As a reproduction scheme 
we employ the Wright-Fisher model (IFlSHER 1930| IWright 193ip . the prototypical model 



of discrete, non-overlapping generations. Since our main concern is how recombination affects 
the speed of adaptation, we assume that all mutations are beneficial. This naturally leads 
us to study evolution in the framework of the infinite-sites model (IKlMURA 19690 : otherwise 



back mutations of beneficial mutations, which are deleterious by definition, should appear 
with nonzero probability. Furthermore, we assume no epistasis among mutations, which will 
be refiected by the multiplicative fitness assignement (see below). As a minimal model with 
the above properties, we study an evolving population with only two loci under selection. 
Each locus is assumed to have infinitely many sites. We assume an initially homogenous 
population and the fitness of the initial genome, or wild-type, is set to unity. 

In line with the assumption of multiplicative fitness effects, the fitness of an individual 
that has rij mutations at locus i compared to the wild-type is exp(r;,iSi + ^2^2). Without loss 
of generality we take S2 > Si. Note that two genotypes with the same number of mutations at 
each locus are not necessarily the same though both have the same fitness exp(sini + S2?t-2)- 
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Since we are only interested in how fast mean fitness increases and not in tfie genealogy, all 
genotypes with the same number of mutations at each locus will be treated as if they were 
the same. 

The population evolves in the following way. Let /t(ni,n2) denote the frequency of all 
genotypes with rii (712) mutations at the first (second) locus in generation t (t = 0, 1, 2, . . .). 
Initially, the population is homogeneous with ft{ni,n2) — (^ni,o<^n2,0) where 5x,y is the Kro- 
necker delta symbol which is 1 (0) ii x = y {x y). By selection, the frequency at generation 
t + 1 on average will change to be 

gSini+S2n2 

//(ni,n2) = z Mni,n2), (1) 

where 

wt= J]e^i"i+^^"Vt(ni,n2) (2) 

ni,n2 

is the average fitness of the population at generation t. 

Mutation can also change the frequency of genotypes. Let go{mi,'m2) denote the proba- 
bility that mutation changes rui sites of the i-th locus (i = 1, 2) of an offspring (mi, m2 > 0). 
Here we implicitly assume that the mutation probability is not affected by the genetic back- 
ground. Note that the probability of having no mutation is included in qq for convenience. 
The frequency change due to both selection and mutation is 

/f(m,^2) = ^ 9o{ni-mi,n2-m2)ft{mi,m2) 

mi,m2>0 

= go{ni - mi,n2 - m2) z /t(mi,m2), (3) 

— Wt 

mi,m2>0 

where go{x,y) with at least one negative argument should be understood to be 0. We 
further assume that mutation does not have any preference for a certain locus, that is, 
go{fni,m2) = go{rn2,mi) for any pair of mi and m2. 

After selection and mutation, two randomly chosen parents mate and beget an offspring. 
Let R{ni,n2\ki, k2;li, I2) denote the probability that the resulting progeny of two individuals 
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with respective genotypes {ki,k2) and {hyh) has the genotype (ni,n2). To be specific, we 
set (0 < r < 1) 

1 — r r 

R{ni, n2\ki, k2] h, h) = ^ (5nifci5?i2fc2 + ^nxh^n2l2) + 2 (^"i'fci^n2'2 + ^n^h^n2k2) ? (4) 

which means that with probabihty 1 — r the two loci of the offspring in question are inherited 
solely from a single parent and with probability r the offspring inherits one locus from one 
parent and the other from the other parent. Note that the case with r = | is generally 
referred to as the free recombination, while r = 1 corresponds to obligate sexuals. 

Since the probability that the randomly chosen parents have genotypes (fci, ^2) and (/i, I2) 
is ft{ki,k2)ft{h,h), the mean frequency after selection, mutation, and recombination is 

f[{ni,n2) = ^ R{ni,n2\ki,k2;lij2)ft{ki,k2)ftih,l2) 

ki,k2,h,l2 

= (l-r)/,>i,n2) + r/«M/f M (5) 

where 

fi'\n,) = J2ftinun2), fi'\n2) = J2ftinun2), (6) 

are marginal frequency distributions of genotypes after the selection and mutation steps with 
ni mutations at locus 1 and n2 mutations at locus 2, respectively. 

Finally, the actual population distribution at generation t + 1 is determined by multino- 
mial sampling using f[{ni,n2) in Equation [5] with the restriction that the population size 



is A^. For simulations, we employ the algorithm explained by Park and Krug (2007) [see 



also Park et al. (2010) for simulations of extremely large populations]. 

The speed of adaptation, or shortly speed, is defined as the rate of increase of the log 
mean fitness, 

. ^ hm ii^^, (7) 

t— >oo t 

where (. . .) denotes an average over independent realizations of evolution with the same 
parameters. In the following, we mainly focus on the dependence of speed on parameters 
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such as the mutation probabihty per generation, the selection coefficient of a single mutation, 
and the recombination probabihty r. 



RESULTS 

Infinite populations: Ahhough the infinite population limit cannot be reached in real 
biological populations for the model we are studying ( IPark et al. 2010"]) . it does provide some 



insight into the adaptation dynamics of finite populations. Furthermore, the deterministic 
nature of the infinite population dynamics renders an analytic approach feasible. 

When the population size is infinite, ^^2) is equal to f[{ni,n2) by the law of large 

numbers. In APPENDIX A, we derive the mean fitness from the recursion relation for the 
frequencies in two special cases, r = and r = 1. From Equations IA7I and lAlll we deduce 
the speed of adaptation as 

- ( n AT ^ r lnGo(e^^*,e^^*) 
Va = v{r = (J, A = 00) = lim 



_ ( , .r , y lnGo(e-^*)+lnGo(e-^*) 
Vs = v{r = 1,N = 00) = lim , (9) 

i— >c»o t 

where subscripts a and s stand for asexuals and obligate sexuals, respectively. Go is the 
generating function for go, defined as 

Go(zi, Z2) = J2 4'4'9oih, h), (10) 

fcl,fc2 

and Go{z) = Go{z,l) = Go{l,z) is the generating function for the marginal distribution 
^o(^) = J2m9oif^^^) = J2m9oi^^ [see also Equation IA5j . Note that we are using the 
symmetry go{ki, /C2) = (70(^2? ^1) introduced earlier, but the generalization to asymmetric 5^0 
is straightforward. 

Since the arguments of Gq in Equation [8] and of Go in Equation [9] increase exponentially, 
the speed is fully determined by the largest possible fitness effect due to a single mutation 
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event. Thus, 

Va = Max{raiSi + n2S2|5'o(^i,'^2) 7^ 0}, (11) 
Vs = M(si+S2), (12) 

where M is the largest possible number of mutated sites at one locus by a single mutation 
event, 

M = Max{n|3n' such that gQ{n,n') = go{n',n) ^ 0}. (13) 

Since, by definition, M is the maximum of all possible n\ and n2 with gijix^ri'i) ^ 0, Vg 
cannot be smaller than Va- Thus, sex is at least not detrimental, though it may have no 
effect depending on the form of qq. For example, if single mutations occur with probability 
U and double mutations involving both loci with probability ?7^, corresponding to 

gQ{k, m) = (l-U - f/^)4+m,o + ^h+m,i + t^^4,i^m,i' (14) 

then Vs = Va = si + S2. On the other hand, if double mutations are forbidden and 

go{k, m) = (1 - U)5k+m,o + ^4+m,i, (15) 

we have Vg = Si + S2 > Va = S2 (recall that we assume S2 > Si). Hence the effect of sex 
significantly depends on the form of go in the infinite population limit. If S2 > Si (strict 
inequality) and if go is as in Equation [T51 beneficial mutations occurring on locus 1 do 
not contribute to the speed of an infinite asexual population. This can be understood in 
the framework of clonal interference as the 'wasting' of weaker beneficial mutations by the 
competition with stronger mutations. On this account, we will later assume Si = S2 when 
we study the adaptation dynamics of finite populations in order to maximize the effect of 
both loci in the asexual case. Note that for si = S2 the mutation scheme in Equation [T5] 
implies Vg = 2va, a two- fold advantage of sex. 

Although we only calculated the speed exactly for the cases r = and r = 1, we now argue 
that the asymptotic speed does not depend on r provided r > 0. Let ii(t) (hit)) denote the 
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maximum number of mutations at locus 1 (locus 2) accumulated up to generation t: 



This definition can be used for finite populations as well, and is closely related to the lead 
of the fitness distribution considered in the traveling wave approach to asexual adaptation 
(IDesai and Fisher 2007t IPark et al. 2010t IFisher 20T3l) . Since the population size is 
infinite, the (marginal) frequency of genomes with £j(t) + M mutations at the ith locus at 
generation t + 1, where M was defined in Equation [131 is nonzero, though it can be extremely 
small. Thus, by recombination, a nonzero fraction of individuals will have 2M + £i(t) +£2{t) 
mutations in total. Combined with the initial condition ^i(O) = ^2(0) = 0, the detectable 
largest fitness at generation t is then M(si + S2)t for all t > 0. Now assume that the speed 
f (r, = 00) for < r < 1 is strictly smaller than M(si + S2). Then, with time t, the ratio of 
the detectable largest fitness to the mean fitness increases as exp((M(si + S2) — v)t). Thus, 
at some point, ft{Mt, Mt) jumps abruptly up to a significant number in a single generation, 
which eventually leads to an abrupt increase of Wt up to the order of exp(M(si + ^2)^). 
Hence, in the long run, the speed becomes M(si + S2) for any r > 0. 

In the above discussion, we argued that the speed does not depend on r once r is nonzero. 
On the other hand, if r is very small, the whole population behaves almost like an asexual 
population for quite some time. Hence, the abrupt jump of fitness mentioned above should 
be observable. To see this phenomenon, we studied the deterministic evolution numerically, 
using the mutation scheme of Equation [T5] with U = 0.1 and si = S2 = 0.02. In Fig. [H 
we show how the mean fitness behaves with time for r = 0, r = 10~^, and r = 1. Even for 
the minute recombination rate of r = 10^^, the mean fitness closely follows the r = 1 curve, 
however with some oscillations. To eludicate the origin of this behavior we need to consider 
how the frequency distribution changes with time. 



£1 (t) = max < n 




= max 




[Figure 1 about here.] 
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In the asexual case (r = 0) the frequency distribution over the number of mutations is 
well described by a Gaussian ( IPark et al. 2010|) . Furthermore, the frequency distribution 



of the obligately sexual population with r = 1 should also be well described by a Gaussian, 
because the generating function is just the product of two generating functions of asexual 
evolution (see Equation lASp . However, for < r ^ 1, the Gaussian may not be a good 
approximation. In Fig. [2] we depict the time evolution of the frequency distribution for 
r = 10^^. Clearly the frequency distribution cannot be approximated by a Gaussian traveling 
wave. Moreover, the shape of the distribution changes with time, which implies that there 
is no time-independent steady state. Rather, the distribution behaves like a 'breathing 
traveling wave' in that the behavior seen in Fig. [2] repeats periodically. The time when two 
peaks become comparable in Fig. [2] corresponds to the abrupt jump of mean fitness alluded 
to above. Further mathematical analysis of this phenomenon seems interesting, but we will 
not pursue it here because it is hardly observable in real, finite populations. 

[Figure 2 about here.] 

Finite populations: We mentioned above that many analytic approaches have been devel- 
oped to find an expression for the speed of adaptation in large asexual populations (IRouziNE et al. 2003 



IDesai and Fisher 200T1 IBrunet et al. 20081 IRouzine et al. 2008p . [Park et al. (2010p 



summarized these developments and compared simulation results with the proposed analytic 



expressions. The approximation of RouziNE et al. (2008 ) turned out to be quite accurate in 



a wide range of parameters. The only disadvantage of this approach is that the speed is ob- 
tained as an implicit function of N (see below). In this section, we will find a mathematical 
formula for the speed of adaptation in sexual populations, using both the suggested formula 
for asexuals and the results for the infinite population dynamics in the previous section. 

For an infinite population, as we have seen above, the detailed form of the mutational 
probability distribution go{ki,k2) affects the speed. However, for plausible values of the 
mutation rate and the selection coefficient such infinite population effects become observable 
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only for unrealistically large populations. To be concrete, it was shown by Park et al. (2010) 
that population sizes of order ~ j^gisoo ^^^^ required to reach the infinite population regime 
for s = 0.02 and U = 10^^. Hence, for a finite population of realistic size, the exact form of 
the mutation probability is not important provided the probability of having n mutations in 
one generation decreases as f/" with U <^ 1. In the following we therefore use Equation [T5] 
and set Si = S2 = s for the reasons mentioned previously. 

We begin with a discussion of the speed for asexual populations. As was illustrated by 



Park et al. (2010), the speed for the asexual version of our model (r = 0) is well approxi- 



mated by the implicit equation 



RBW / RBW \ / 3TT 

IniV - ^ In^ ^ + 1 - In J-WT^n-T-WWTTTT-TT, (17) 



2s2 V eUs J y t;^BWin(yRBW/(^5)) 
where the subscript a in Va refers to the asexual population, the superscript RBW refers to 



the authors of RouziNE et al. (2008), and e ^ 2.718182 is the base of the natural logarithm. 

Since the approximation of the fitness distribution by a continuous traveling wave was 
used to derive Equation [T71 it should not be surprising that the discrepancy between the- 
ory and simulation becomes relatively large when the size of population is small enough to 
realize the strong-selection weak-mutation (SSWM) regime, where the population is mostly 
monomorphic. Based on this observation, there is room for improvement of the approxima- 
tion in an ad-hoc way as follows: First we note that the first (second) term in Equation [17] 
is dominant when the speed is high (low). Thus, when the population size is small, we can 
neglect the first term. Recalling that the speed in the SSWM regime is given by Va = 2NU 



( IPark et al. 2010jl . we modify Equation [IT] as 

IniV^^fln^^ + l) +ln^ (18) 
2s2 V eUs J 2s^U ^ ' 

which keeps the large speed behavior unchanged and enforces the SSWM result for small 
speeds. In Fig. [3] we show that Equation [18] provides a more accurate approximation to the 
speed obtained from simulations with U = 10~^ and s = 0.01 than Equation [T71 
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[Figure 3 about here.] 

Now we move on to the speed of sexual populations. At first, let us start from the case of 
r = 1 whose infinite population limit allows for an exact solution. As we show in APPENDIX 
A, the evolutionary dynamics of an infinite population with r = 1 can be viewed as the 
independent evolution of each locus with the marginal mutation probability go{k). That 
is, we can divide the evolutionary dynamics into two independent asexual populations with 
reduced mutation probability and the speed of the sexual population is obtained by simply 
adding the speeds of these two virtual asexual populations. Within the mutation scheme 
given by Equation [15] with selection coefficients si = S2 = s this implies that 

V, (r = l,U) = 2va iU/2) (19) 

for sufficiently large populations. Interestingly, Equation [19] is trivially valid in the SSWM 
regime where the speed is linear in U and v{r,U) ~ 2NUs^ irrespective of r. Since Equa- 
tion [19] accurately estimates the speed for very small and very large populations, it is likely 
that Equation [19] is a good approximation for any population size. Indeed, as we show in 
Fig. m Vs{r = 1,U) is well approximated by twice Va{U/2) for any population size. The 
parameters we have used in these simulations are U = 10^^ and s = 0.01. With the help of 
Equation [TS] we may thus approximate the speed Vg of the obligately sexual population as 

In AT ^ ^ f In^ ^ + iK In (20) 
4s2 V eUs J 2s^U ^ ' 

[Figure 4 about here.] 

It is clear that for < r ^ 1 there should be a regime where Equation [19] cannot 
approximate the speed accurately. To see this deviation, we simulated populations with 
various r (Fig. Hj). It turns out that Equation [19] is still a good approximation for r > 10^^. 
In particular, the speed for r = 0.1 is hardly discernible from that for r = 1 for all population 
sizes. The deviation starts to be significant for r = 10~^. For comparison, we also plot 
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Vs{r = 0, U)/va{U/2) or equivalently Va{U)/va{U/2) in Fig. HI which should approach to 1 in 
the infinite population limit (IPark et al. 2010|) . For > 10^, where NU \a{N s) becomes 
larger than 1, Vs{r,U)/va{U/2) starts to increase though very slowly and Vs{r,U) becomes 
significantly larger than Vs{r = 0,U) = VaiU). Note that for asexual populations clonal 
interference sets in around NU\n{Ns) ~ 1 (IWlLKE 20041 IPark et al. 2010p . That is, as 
soon as clonal interference becomes relevant, even a small amount of recombination leads to 
a significant speedup of adaptation, in agreement with the FM mechanism. The fact that 
the ratio Vg/va continues to rise monotonically with for all cases with r > in Fig. H] is 
consistent with the two-fold advantage predicted by the infinite population analysis. 

When measuring the speed of adaptation in our simulations, a useful consistency check 
was provided by the Guess relation 



/ 1 

\ i=i 



X. ) , (21) 



Stat 



where Xi is the relative fitness of i-th individual in the infinite time limit, 

X. = hm (22) 

t^oo w[t) 

and (■)stat signifies an average over the stationary measure of Xi- Equation [21] was orig- 
inally established for asexual populations undergoing discrete generation (WF) dynam- 
ics fIGUESS 1974at IGUESS 1974bll . In APPENDIX B, we prove that the relation holds for 



sexuals as well, and in Fig. [5l we numerically confirm its validity. The two terms on the 
right hand side of Equation [21] represent the increase in population fitness due to mutation 
and selection, respectively. Recombination affects the speed of adaptation only indirectly 
through its effect on the relative fitnesses Xi- Note that the Guess relation should hold even 
if one uses discrete-time, overlapping generation models such as the Moran model. 

[Figure 5 about here.] 

Finally, we analyze the difference in the number of beneficial mutations acquired by the 
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two loci. We quantify this difference as 



a = lim 



(23) 



t 



where ii and £2 are defined in Equation [161 We will refer to a as the mutation number 
imbalance (MNI). To discern the MNI of asexuals from that of sexuals, we will add subscripts 
a and s, for asexual and sexual populations, respectively. In infinite populations each locus 
accumulates the same number of mutations, hence this study is meaningful only for finite 
populations. 

For asexual populations, an approximation for a can be obtained by comparing the 
origination processes at the two loci, which count the mutations that are present in some 
individuals of the population at time t and that are destined to eventually go to fixation 
( IGlLLESPlE 1993| IGlLLESPlE 19941 IPark and Krug 20071) . Denoting the number of such 
mutations at locus i by ki{t), we assume that 1) the difference between ki{t) and the lead 
ii{t) > ki{t) remains bounded in the long time limit, 2) the total number of mutations 
Mp = ki + k2 in the origination process increases at the same rate as the mean number of 
mutations, Mp{t) ^ Vat/s for large t, and 3) each new mutation appearing in the origination 
process chooses one of two loci with equal probability. Assumptions 1) and 2) reflect the 
existence of a steady state and have been verified in simulations (IPark and Krug 2007p . 
and assumption 3) is a consequence of the symmetry between the two loci. By assumption 
3), the probability that there are mi mutations at locus 1 and m2 = Mp{t) — mi mutations 
at locus 2 is given by 



Since the mean of mi is (mi) = Mp/2 and its variance is ((mi — (mi))^) = Mp/4, we can 
calculate aa, invoking the assumption 1), as 



In Figure ini we compare a a to Va/s for U = 10 ^ and s = 0.01, which shows an excellent 




(24) 




(25) 



agreement. 
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[Figure 6 about here.] 

Recombination changes the behavior of the MNI substantially. As can be seen in Figure O 
once the population size is in the regime of clonal interference the MNI decreases abruptly, 
then remains almost constant for a wide range of population sizes. Even a small amount 
of recombination efficiently equalizes any major fitness difference between the two loci by 
creating competitively superior recombinants in which both loci have high fitness. 

DISCUSSION 

The Fisher-Muller mechanism for the evolutionary advantage of sex is based on the slowing 
down of asexual adaptation due to clonal interference, which is alleviated by the recombi- 
nation of high fitness genotypes. While much recent theoretical work has been devoted to 
quantifying the speed of adaptation in asexuals, the speedup that can be achieved through 
recombination has been explicitly addressed only in a few studies (see below). In the present 
article we take a step in this direction by providing a detailed analysis of a simple, yet bio- 
logically meaningful model in which recombination occurs between two loci, each of which 
can harbor an unlimited number of linked beneficial mutations. Our analysis shows that the 
advantage of sex becomes significant in the parameter regime where clonal interference plays 
an important role in asexual populations. In our two-locus model, the adaptation speed of 
sexual populations is about twice as large as that of the corresponding asexual populations 
for a wide range of recombination rate, an effect that we refer to as a two-fold advantage of 
sex. In the remainder of this section we discuss the extension of our model to L loci and 
compare our findings to related previous work. 

More than two loci: It is natural to surmise that the factor of two arises in our model 
because we are considering two loci, and that the speed increase should in general by pro- 
portional to the number of loci. Indeed, this turns out to be true if we use a 'communal' 
recombination scheme where the gene of each locus is collected from the whole population 
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rather than from the two parents, once the genome of the offspring is constructed by recom- 
bination. In a three-locus model with this recombination scheme, the frequency distribution 
of next generation is sampled from 



where f^ irii) is the marginal frequency distribution for having mutations at locus i after 
the (deterministic) selection and mutation steps (compare to MODELS). It is a straight- 
forward extension of the calculation in APPENDIX A to show that the infinite population 
dynamics for r = 1 is again divided into three independent evolutions of each locus with 
marginal mutation probabilities just as in the two-locus case. In general, if we consider a 
model system with L loci within the communal recombination scheme mentioned above, the 
evolution is a superposition of L independent evolutions at each locus, and there is an L-fold 
advantage of sex in the infinite population. 

To see if this L-fold advantage persists for finite populations, we performed simulations 
of the three-locus model. As in the two-locus case, we expect that 



for sufficiently large r. Indeed, we observe that the simulations are consistent with Equa- 
tion [27] for a wide range of parameter values. In Fig. [71 we depict Vs{r,U)/va{U/3) as a 
function of for U = 1.5 x 10^^ and s = 0.01 with varying r. As in the two-locus model, 
the advantage of sex becomes significant when clonal interference is important in the corre- 
sponding asexual populations. 



We also studied the mutation number inbalance in the three-locus model. We slightly modify 
the definition of the MNI as the difference between the maximum and minimum numbers 
of accumulated mutations at all loci, which reduces to the definition of Equation [23] for the 



f[{ni,n2,ns) = (1 - r)//'(ni, rig, ng) + r/f ^^2)/^ Vs), 



(26) 




(27) 



[Figure 7 about here.] 
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two-locus model. In Fig. |8l we depict the MNI for the three- locus model with U = 1.5 x 10^^ 
and s = 0.01 for various r. Like the MNI of the two-locus model, the MNI for the asexuals 
increases with while slowly decreasing for sexuals. Again, the qualitative difference be- 
tween sexuals and asexuals becomes significant in the regime where clonal intereference is 
important. 

[Figure 8 about here.] 

Relation to previous studies: The quantitative analysis of the speed of evolution of sexual 
populations compared to that of asexual populations, when both evolve on the same non- 
epistatic fitness landscape with the same beneficial mutation rate per genome, has a long his- 
tory (ICrqw and KiMURA 1965t IMaynard Smith 1968j ICrqw and Kimura 1969t IMaynard Smith l£ 
IFelsenstein 1974t IMaynard Smith 19781 IKim and Orr 20051) . When the number of 
accessible beneficial mutations is finite, the relevant quantity is the time for all beneficial 



mutations to be fixed. In this context, Maynard Smith (1971) analyzed the fixation time 



for sexual and asexual populations evolving on a fitness landscape with / loci under selection. 
Each locus has two alleles, one of which confers a beneficial fitness effect in a non-epistatic 
fashion. For sexual populations, the linkage among loci was assumed weak. Using a rough 



approximation, Maynard Smith (1971) argued that the time for completing evolutionary 
changes in asexual populations is / times longer than that in sexual populations for suffi- 
ciently large population size, which implies an /-fold advantage of sex similar to what we 



found in our study. While the analysis by Maynard Smith (1971) is fairly crude and (as 



conceded by the author) actually not consistent with the simulation results presented in the 
same paper, the conclusion that the advantage of sex becomes stronger with an increasing 
number of loci under selection is in qualitative agreement with our results, as well as with 



the related work of KiM and Orr (2005). 



Recent studies of the speed of sexual populations in the context of the FM mechanism 
have mostly focused on the case with an infinite supply of beneficial mutations ( INeher et al. 2010} 
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IRouziNE and Coffin 2010} IWeissman and Barton 201"2|) . exploiting the mathematical 

progress in treating the spreading of beneficial mutations as a Gaussian traveling wave ( ITsiMRlNG et al. 199 

IRouziNE et al. 20031 IDesai and Fisher 20071IRouzine et al. 20081 IPark et al. 20101 [Good et al. 2012 



IFlSHER 20T3l) . [RouziNE and Coffin (2010[) studied how recombination speeds up adapta- 



tion when there is standing variation of beneficial mutations. Neher et al. (2010) studied 



the speed of adaptation of large facultatively sexual populations, starting from a monomor- 



phic state. Similar to our results, Neher et al. (2010) found a regime of intermediate re 



combination rates where the speed increases logarithmically with population size, however 



with a prefactor that varies quadratically with r. Although RouziNE and Coffin (2010) 



and Neher et al. (2010) investigated the adaptation dynamics of sexual populations with 
an (effectively) infinite supply of beneficial mutations, their results cannot be directly com- 



pared to ours. This is because the model genomes of RouziNE and Coffin (2010) and 



Neher et al. (2010) assume no or weak linkage between beneficial mutations, whereas in 



our model mutations in the same locus are tightly linked. A related study with an explicit 



genetic map was recently presented by Weissman and Barton (2012). 

The two-locus genome considered in this paper can be viewed as a simple example of a 
modular genomic architecture, where recombination occurs between modules but not within 



a modul. Watson et al. (2011) have pointed out that such a modular structure induces 



a strong benefit for sexual reproduction when there is sign epistasis within the modules 
and different modules contribute independently to fitness. It would therefore be of inter- 
est to extend our approach to include a tunable degree of epistasis interactions within the 



loci. Following Watson et al. (2011), such interactions should affect not only the speed of 



adaptation but also the set of genotypes that can be reached at all by the population. 
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APPENDIX A: INFINITE POPULATION DYNAMICS FOR ASEXUALS (r = 0) AND 

OBLIGATE SEXUALS (r = 1) 

Let Ft{zi,Z2) denote the generating function for the frequency distribution at generation t, 
which is defined as 

Ft{zr,Z2)^ J2^i'^TMnun2). (Al) 

Since the fitness landscape is multiplicative, the mean fitness at generation t can be found 
from Ft through 

Wt = Ft{e'\e''). (A2) 

Likewise, we introduce the generating function for // in Equation [H which is obtained from 
Ft according to 

^ Y: -TzTmnun,) = ''yi^:;'^;;"^^ (A3) 

711 ,712 t\, 1 J 

Since fj^ in Equation [3] is the convolution of go and //, the generating function for ft is the 
product of F/ and Go{zi,Z2) (for the definition of the latter quantity see Equation [TOj) . 

Using that ft+i is the same as fj for infinite populations, we obtain an iterative evolution 
equation for Ft that reads 

Ft+M,Z2) = (l-r)Ff(zi,Z2)+rF,'^(l,^2)i^r(^i,l) (A4) 

n , Ftiz,e'\z2e^') , ~ ^Kzle^e^~ ^Ke%^2e^ 
= (1 — rjGoizi, zo) -, ^ h rGn Zi z ^LrnZo ; 

where Fj^^zi, Z2) = Go{zi, Z2)Ff{zi, Z2) is the generating function of and 

Go{z) = J2 ^''9o{ki, k2) = J2 ^''9oih, k2). (A5) 

fci,fc2 A;i,fc2 
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For r = 0, Equation lA4l can be solved by iterating backwards until t = 0, that is (IPark and Krug 2007P 

Ft{zi, Z2) = Go{zi, Z2)—- — — — -— = Go{zi, Z2) 



^ F,{z^e'^\z2e'^-) yr G^{z^e^^^^Z2e^ 

with ^0(^1,^2) = 1 for the homogeneous initial condition. Thus the mean fitness at genera- 
tion t is 

«7t = Go(e^^*,e^^*). (AT) 

When r = 1, Equation IA4I suggests that the generating function is the product of two 
functions such that 

F,iz^,Z2)=Fl{z^)F^iz2). (A8) 

With the above ansatz, we can find an evolution equation for F^{z) = 1 or 2) from 
Equation [M 

FU^) = Go{z)^-^, (A9) 

which is exactly the evolution equation for an asexual population with mutation generating 
function Go{z). The solution of Equation IA9I is 

fi G„(e-) 

where we have used the initial condition ^0(2:1, 2:2) = -^0(2^) = 1- One can easily check that 
Equation IA8I with Ft{z) given by Equation IA9I actually solves Equation IA4I for r = 1 by 
substitution. Hence the mean fitness at generation t for r = 1 is 

= Ft{e'\e'^) = Go(e^^*)Go(e^^*). (All) 

APPENDIX B: GUESS RELATION IN THE PRESENCE OF RECOMBINATION 
In this APPENDIX, we will show that for evolution on multiplicative, non-epistatic fitness 
landscapes the Guess relation (Equation [2Ti) is valid even in the presence of recombination. 
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Let Wi(t) be the fitness of the i-th individual at generation t, w{t) the mean fitness 
of the population, w{t) = '^i'Wiit)/N, and Xi{t) the relative fitness of i-th individual, 
Xi{t) — Wi{t)/w{t). We will assume that Xi{t) approaches a well-defined steady state as t 
goes to infinity. We take each individual to be characterized by a genome with L loci, each 
of which has inifitely many sites. The contribution of a locus to fitness is denoted by Zn 
{n — 1, . . . , L) and the fitness of an individual with such a genome is w — Ylt=i ^n- In the 
following, Zn will be called CFn, meaning the Contribution to Fitness of the nth-locus. If a 
mutation hits the n-th locus, CFn changes from Zn to z'^ = ZnVn, where Vn is drawn from a 
given probability distribution that may vary from locus to locus but does not depend on Zn 
or the generation. If Vn is larger (smaller) than 1, the mutation is beneficial (deleterious). 
In this APPENDIX, the explicit form of the probability distribution for Vn does not need to 
be specified. 

We will use the vector notation z= {zi, zl) for fitness vectors with L elements. Assume 
that there are N individuals and the CFn of individual i is Zi^n- The corresponding fitness 
vector is denoted by Zi. At first, we will calculate the expected mean log- fitness at generation 
t + 1 assuming that Xi{t) and Zi{t) are given. 

By selection, the probability density that the fitness vector of an offspring is z is 



where S{x) is the L-dimensional Dirac delta function. Let g{v) be the probability density 
that a mutation event changes the CFn of an offspring by Vn {zn z'^ — ZnVn) for all n's. 
Then due to mutation, the expected frequency becomes 



where z' ®v denotes the vector with elements z'^Vn- 

Next we consider recombination. Let R{z\zi, Z2) be the probability density that offspring 
resulting from the recombination of two parents with fitness vectors zi and Z2 has fitness z. 




(Bl) 



1=1 




(B2) 
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In general, we can write R in the form 

L 

R{z\Zi, Z2) = ^piS) JJ (5 (Z„ - Zs(n),n) , (B3) 
5 n=l 

where S runs over all possible outcomes of recombination and p{S) is the probability of this 
event. Here S{n) = 1 (2) if locus n is inherited from parent 1 (2), hence the total number of 
possible outcomes is 2^. Then the final probability density becomes 

/(^ = y dZidZ2R{z\zi,Z2)fm{zi)fm{z2). (B4) 

Now we will calculate the expected log-fitness of a randomly chosen individual at gener- 
ation if: + 1 for given Wi{t). Since the probability density that a randomly chosen individual 
at generation t + 1 has fitness vector z is f{z) given in Equation IB4I and the corresponding 
fitness is w = Yin ^n, the quantity we want to calculate is 

/ = y (iwlntu prob(w) = J cizln ^JJ^ 2;^^ /(i*), (B5) 
where prob(w) is the probability density that an individual has fitness w at generation t + 1 



for given Wi{t). GuESS (1974b) showed that as long as there is a well-defined steady state 



the speed v can be calculated as 



ln^^^^\ = /ln^') = (/)-(lnW.), (B6) 



w{t) I \ W2, 

where Wi (W2) is the fitness of a randomly chosen individual at generation t + 1 (t) and (. . .) 
signifies the average over the steady state distribution. Loosely speaking, the above relation 
can be understood as follows: Since the dynamics is symmetric under permutations of the 
population index, the steady state must have this permutation symmetry as well. Hence the 
expected log-mean fitness at steady state should be the same as the expected log-fitness of a 
randomly chosen individual. Moreover, at stationarity the speed can be calculated from the 
difference of log-fitness between two consecutive generations. Thus Equation IB6I follows. 
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From Equations EH |B2l [B3l and EH we get 



I = [ dz\n(l[zJJ2 

\n=l J i,j 



'J 



dvidvjR{z\zi ® Vi, Zj ® Vj)g{vi)g{vj) 
^'^^^^^^^^ I dvidvjg{vi)g{vj) f dzln | J]^ ^^n | R{z\zi ® Vi, zj (g) vj) 

i,j ^ \n=l J 



i,j S n=l 

^PiS) ^ I ln(2;s(„),„) + j dvidvjg{vi)g{vj) \n{vs(„),n) j , (B7) 



S n=l \ i,j 

where S{n) in the subscript of z and v can be either i or j. Since 



dvidvjg{vi)g{vj)\nvi^n = / dvidvjg{vi)g{vj)\nvj^n = / (it75(({;) lnf„, (B9) 



the summation and integral in the parentheses of the last line of Equation IB7I do not depend 
on S. Due to the normalization YlsPi^) ~ S^^ 

^ = ^^ln(X.(t)«;(t)) + (InV) = ^ ln(X,(t)) + (In^) + ln«;(t), (BIO) 

i i 

where V is the total effect of fitness increase by a single mutation event, V = Yln'^n, and 
rin-^*." ~ "was used. Hence the speed becomes 

where Xi is the relative fitness of the individual i at steady state. Note that the formula 
does not depend on the explicit form of the recombination operator. If we use Equation [15] 
for the mutation scheme, (In V) = Us. 

The application of the above procedure to the Moran model is straightforward. Hence, 
the Guess relation is valid for discrete time models regardless of recombination, once the 
fitness landscape is multiplicative. 
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List of Figures 

Mean fitness Wt of tlie infinite population model as a function of time for 
r = 0, r = 10~^, and r = 1 (from bottom to top) with f/ = 0.1 and s = 0.02. 
As argued in the text, the speed does not depend on r once r is nonzero. . . 
Frequency distribution of the total number of mutations for the infinite pop- 
ulation model at generations 895, 900, 905, . . . , 930 (left to right) with param- 
eters r = 10-^ U = 0.1, and s = 0.02 

Speed of adaptation of finite asexual populations as a function of on a 
double logarithmic scale for U = 10~^ and s = 0.01. The numerical solutions 
of Equation [T7] and Equation [18] are drawn for comparison with the simulation 
data. As anticipated, the ad-hoc modification (Equation [T8l) provides a more 

accurate estimate 

Ratio of the sexual adaptation speed, Vs{r, U), to the asexual speed at half mu- 
tation rate, Va{U/2), as a function of population size A^, where Vs is obtained 
from simulations and Va is computed using Equation [181 Recombination rates 
are r = (empty reverse triangle), 10~^ (empty square), 10~^ (filled trian- 
gle), 10~^ (empty triangle), 10~^ (filled circle), 10~^ (empty circle), and 1 (filled 
square) from bottom to top, and U = 10~^ and s = 0.01 are used throughout. 
The scaling relation in Equation [T9] predicts that Vs{r, U)/va{U/2) = 2. Note 
that two datasets for r = 0.1 (empty circle) and r = 1 (filled square) are 

indiscernible 

Numerical verification of the Guess relation. The figure compares the loga- 
rithmic mean fitness divided by generation time, {lnw{t)) /t, to the right hand 
side of Equation [211 The data are obtained from simulations with r = 10~^, 
U = 10~^ and s = 0.01. For t — )• oo, both curves intersect the ordinate at the 

same point, which equals the asymptotic speed of adaptation 

The mutation number imbalance (MNI) a vs. population size A^ for r = 
(empty reverse triangle), 10~^ (empty square), 10"'^ (filled triangle), 10^^ 
(empty triangle), 10~^ (filled circle), 10~^ (empty circle), and 1 (filled square) 
from top to bottom. Other parameter values are U = 10~® and s = 0.01. 
The symbols for r > 10~^ essentially overlap. For comparison, the analytic 

prediction aa = Va/s for asexuals is drawn as a line 

Ratio of the sexual adaptation speed in the three- locus model, Vs{r, U), to the 
asexual speed Va at mutation rate U/3 as a function of population size A^, 
where Vs is obtained from simulations and Va is computed using Equation [181 
Recombination rates are r = (empty reverse triangle), 10~^ (empty square), 
10~^ (filled triangle), 10~^ (empty triangle), 10~^ (filled circle), 10~^ (empty 
circle), and 1 (filled square) from bottom to top, and U = 10~^ and s = 
0.01 are used throughout. The scaling relation in Equation [27] predicts that 
Vs{r,U)/va{U/3) = 3. Note that two datasets for r = 0.1 (empty circle) and 
r = 1 (filled square) are indiscernible 
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Mutation number imbalance (MNI) for the three-locus model with recombi- 
nation rates r = (empty reverse triangle), 10~^ (empty square), 10~^ (filled 
triangle), 10~^ (empty triangle), 10~^ (filled circle), 10~^ (empty circle), and 1 
(filled square) from top to bottom. Other parameter values are f/ = 1.5 x 10~^ 
and s = 0.01. The symbols for r > 10~^ essentially overlap 
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Figure 1: Mean fitness Wt of tlie infinite population model as a function of time for r = 0, 
r = 10~^, and r = 1 (from bottom to top) with U = 0.1 and s = 0.02. As argued in the 
text, the speed does not depend on r once r is nonzero. 
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Figure 2: Frequency distribution of the total number of mutations for the infinite population 
model at generations 895, 900, 905, . . . , 930 (left to right) with parameters r = 10^^, U = 0.1, 
and s = 0.02. 
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Figure 3: Speed of adaptation of finite asexual populations as a function of on a double 
logarithmic scale for U = 10~^ and s = 0.01. The numerical solutions of Equation [T71 and 
Equation [18] are drawn for comparison with the simulation data. As anticipated, the ad-hoc 
modification (Equation [18]) provides a more accurate estimate. 
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Figure 4: Ratio of the sexual adaptation speed, Vs{r, U), to the asexual speed at half mutation 
rate, Va{U/2), as a function of population size N, where Vg is obtained from simulations and 
Va is computed using Equation [181 Recombination rates are r = (empty reverse triangle), 
10"^ (empty square), 10"'^ (filled triangle), 10~^ (empty triangle), 10~^ (filled circle), 10~^ 
(empty circle), and 1 (filled square) from bottom to top, and U = 10~^ and s = 0.01 are 
used throughout. The scaling relation in Equation \W\ predicts that Vs{r,U)/va{U/2) = 2. 
Note that two datasets for r = 0.1 (empty circle) and r = 1 (filled square) are indiscernible. 
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Figure 5: Numerical verification of tlie Guess relation. The figure compares the logarithmic 
mean fitness divided by generation time, (lnw(t))/t, to the right hand side of Equation m\ 
The data are obtained from simulations with r = 10~^, U = 10~^ and s = 0.01. For t — )■ oo, 
both curves intersect the ordinate at the same point, which equals the asymptotic speed of 
adaptation. 
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Figure 6: The mutation number imbalance (MNI) a vs. population size N for r = (empty 
reverse triangle), 10~^ (empty square), 10~^ (filled triangle), 10~^ (empty triangle), 10~^ 
(filled circle), 10~^ (empty circle), and 1 (filled square) from top to bottom. Other parameter 
values are U = 10~^ and s = 0.01. The symbols for r > 10~^ essentially overlap. For 
comparison, the analytic prediction ctq = Va/s for asexuals is drawn as a line. 
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Figure 7: Ratio of the sexual adaptation speed in the three- locus model, Vs{r,U), to the 
asexual speed Va at mutation rate U/3 as a function of population size N, where Vs is 
obtained from simulations and Va is computed using Equation [181 Recombination rates are 
r = (empty reverse triangle), 10~^ (empty square), 10"'^ (filled triangle), 10~^ (empty 
triangle), 10~^ (filled circle), 10~^ (empty circle), and 1 (filled square) from bottom to top, 
and U = 10^® and s = 0.01 are used throughout. The scaling relation in Equation [271 
predicts that Vs{r,U)/va{U/3) = 3. Note that two datasets for r = 0.1 (empty circle) and 
r = 1 (filled square) are indiscernible. 
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Figure 8: Mutation number imbalance (MNI) for the three-locus model with recombination 
rates r = (empty reverse triangle), 10~^ (empty square), 10~^ (filled triangle), 10~^ (empty 
triangle), 10^^ (filled circle), 10~^ (empty circle), and 1 (filled square) from top to bottom. 
Other parameter values are U = 1.5 x 10~^ and s = 0.01. The symbols for r > 10~^ 
essentially overlap. 
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